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Abstract. We study the zero-temperature phase diagram of the Ji -J2 
Heisenberg model for spin-1 particles on an infinite square lattice interacting via 
nearest-neighbour (Ji = 1) and next-nearest-neighbour (J2 > 0) bonds. Both bonds 
have the same XXZ-iype anisotropy in spin space. The effects on the quasiclassical 
Neel-ordered and coUinear stripe-ordered states of varying the anisotropy parameter 
A is investigated using the coupled cluster method carried out to high orders. By 
contrast with the spin-i case studied previously, we predict no intermediate disordered 
phase between the Neel and collinear stripe phases, for any value of the frustration 
J2/J1, for either the z-aligned (A > 1) or a;?;-planar-aligned (0 < A < 1) states. The 
quantum phase transition is determined to be first-order for all values of J2/J1 and 
A. The position of the phase boundary J|(A) is determined accurately. It is observed 
to deviate most from its classical position = 5 (fo^' values of A > 0) at the 
Heisenberg isotropic point (A = 1), where J|(l) = 0.55 ± 0.01. By contrast, at the 
XY isotropic point (A = 0), we find J|(G) = 0.50 ± 0.01. In the Ising limit (A 00) 
J| 0.5 as expected. 



PACS numbers: 75.10.Jm, 75.30.Gw, 75.30.Kz, 75.50.Ee 



1. Introduction 

In a recent paper [1] we have used the coupled cluster method (CCM) [2-4] to study the 
influence of spin anisotropy on the ground-state (gs) magnetic ordering of an anisotropic 
version (viz., the J^^^-J^^^ model) of the well-known J1-J2 model on the infinite 
two-dimensional (2D) square lattice, described below, for particles with spin quantum 
number s = |. In the present paper we further the investigation of the J^^^-J^^^ 
model by replacing the spin-| particles by particles with s = 1. 

The main purpose of the previous paper was to examine carefully the role of spin 
anisotropy in tuning the quantum fluctuations that play such a key role in determining 
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the quantum phase diagram of the pure (spin-isotropic) J1-J2 model that has become 
an archetypal model for discussing the subtle interplay between the effects due to 
quantum fluctuations and frustration, as discussed below. While increasing the spin 
quantum number s is, of course, expected to reduce the effects of quantum fluctuations, 
new and unexpected phenomena may also arise. Thus, a well-known example of such 
new behaviour emerging when s is increased is the appearance of the gapped Haldane 
phase [5] in s = 1 one-dimensional (ID) chains, which is not present in their s — ^ 
counterparts. 

The basic (spin-isotropic) J1-J2 model with nearest-neighbour (NN) and next- 
nearest-neighbour (NNN) antiferromagnetic exchange interactions, of strengths Ji 
and J2 respectively, has been extensively studied both theoretically [6-20] and 
experimentally [21-24]. Many of the earlier studies were motivated, at least in part, 
by the hope of shedding light on the possible link between antiferromagnetism and the 
onset of superconductivity at high temperature in the doped cuprate materials whose 
undoped precursors are seemingly well described by the s — ^ version of the J1-J2 model 
on the square lattice in two dimensions [8, 25-27] . The recent discovery of several other 
quasi-2D materials that are realizations of the J1-J2 model, has only served to extend 
the theoretical interest in the model. 

Some of the actual magnetic compounds that can be well described by the 
s = I J1-J2 model are La2Cu04 [27] for small values of J2/ Ji, and Li2VOSi04 
and Li2VOGe04 [21,22] for large values of J2I J\- Other such materials include 
the compounds VOM0O4 [23] and Pb2VO(P04)2 [24]. The compound VOM0O4 is 
interesting because its exchange couplings appear to be more than an order of magnitude 
larger than those of Li2VOSi04, even though the structures of the two compounds are 
closely related. Similarly, the compound Pb2VO(P04)2 also has a structure closely 
related to that of Li2VOSi04, but it appears to have a ferromagnetic NN exchange 
coupling (Ji < 0) frustrated by an antiferromagnetic NNN exchange coupling (J2 > 0), 
with \ J2l J\\ ~ 1-5. By contrast, although all of the other compounds mentioned above 
are also examples of quasi-2D frustrated spin-| magnets, they have NN and NNN 
exchanges that are both antiferromagnetic. 

For the past few decades, a great deal of attention has also been devoted to magnetic 
materials with spin-1 ions, such as the linear chain systems including CsNiCls [28] 
with a weak axial anisotropy, CsFeBra [29] with a strong planar anisotropy and the 
complex materials NENP (Ni(C2H8N2)2N02(C104)) [30] with a weak planar anisotropy 
and NENC (Ni(C2H8N2)2Ni(CN4)) [31] with a strong planar anisotropy; as well as 
the 2D Heisenberg antifcrromagnet K2NiF4 [32]. The spin gaps observed in CsNiCls 
and NENP are believed to be examples of the integer-spin gap behaviour predicted by 
Haldane [5]; whereas half-odd-integer spin sytems are gapless. Another new spin-gapped 
material is the 2D triangular lattice antifcrromagnet NiGa2S4 [33] which, it has been 
argued [34,35], may be a "spin nematic" [36]. It is clear, therefore, that the theoretical 
study of 2D spin-1 quantum magnets is worthy of pursuit. 

In this context we note the recent discovery of superconductivity with a transition 
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temperature at ~ 26 K in the layered iron-based compound LaOFeAs, when doped 
by partial substitution of the oxygen atoms by fluorine atoms [37], La[Oi_xFx]FcAs, 
with X ~ 0.05-0.11. This has been followed by the rapid discovery of superconductivity 
at even higher values of Tc (> 50 K) in a broad class of similar doped quaternary 
oxypnictide compounds. Enormous interest has thereby been engendered in this class of 
materials. Of particular relevance to the present work are the very recent first-principles 
calculations [38] showing that the undoped parent precursor material LaOFeAs is well 
described by the spin-1 J1-J2 model on the square lattice with Ji > 0, J2 > 0, and 
J2/J1 2. Broadly similar conclusions have also been reached by other authors [39]. 

Many of the above quasi-2D magnetic materials, and many others like them, 
display interesting gs phases, often with subtle quantum phase transitions between them. 
Gcnerically, the interplay between reduced dimensionality, competing interactions and 
strong quantum fluctuations, seems to generate a number of new states of condensed 
matter with orderings that differ from the usual states of quasiclassical long-range order 
(LRO). Thus, for high-temperature superconductivity, for example, Anderson [25] has 
suggested that quantum spin fluctuations and frustration due to doping could lead to 
the collapse of the 2D Neel-ordered antiferromagnetic phase present at zero doping, 
and that this could be a mechanism that drives the superconducting behaviour. This, 
and many similar experimental observations for other magnetic materials of reduced 
dimensionality, has intensified the study of order-disorder quantum phase transitions. 
Thus, low-dimensional quantum antiferromagnets have attracted much recent attention 
as model systems in which strong quantum fiuctuations might be able to destroy 
magnetic LRO in the ground state (GS). In the present paper we consider a system 
of — > 00 spin-1 particles on a spatially isotropic 2D square lattice. 

The isotropic Heisenberg antiferromagnet with only nearest-neighbour (NN) bonds, 
all of equal strength (Ji > 0), exhibits magnetic LRO at zero temperature on such 
bipartite lattices as the square lattice considered here. A key mechanism that can then 
serve to destroy the LRO for such systems (with a given lattice and spins of a given spin 
quantum number s) is the introduction of competing or frustrating bonds on top of the 
NN bonds. The interested reader is referred to [40,41] for a more detailed discussion of 
2D spin systems in general. 

In this context, and as we have already noted above, an archetypal frustrated model 
of the above type that has attracted much theoretical attention in recent years is the 2D 
J1-J2 model on a square lattice with both NN and NNN antiferromagnetic interactions, 
with strength Ji > and J2 > respectively. The NN bonds Ji > promote Neel 
antiferromagnetic order, while the NNN bonds J2 > act to frustrate or compete with 
this order. All such frustrated quantum magnets continue to be of great theoretical 
interest because of the possible spin-liquid and other such novel magnetically disordered 
phases that they can exhibit (and see, e.g., [42-44]). 

The properties of the s = 1/2 J1-J2 model on the 2D square lattice are well 
understood in the fimits when J2 = or Ji = 0. For the case when J2 = 0, and 
the classical GS is perfectly Neel-ordered, the quantum fiuctuations are not sufficiently 
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strong enough to destroy the Neel LRO, although the staggered magnetization is reduced 
to about 61% of its classical value. The opposite limit of large J2 is a classic example [8] 
of the phenomenon of order by disorder [45,46]. Thus, in the case where Ji — with 
J2 7^ and fixed, the two sublattices each order antiferromagnetically at the classical 
level, but in directions which are independent of each other. This degeneracy is lifted 
by quantum fluctuations and the GS becomes magnetically ordered coUinearly as a 
stripe phase consisting of successive alternating rows (or columns) of parallel spins. It 
is by now also widely accepted that the s — 1/2 J1-J2 model exhibits the above two 
quasiclassical antiferromagnetic phases with LRO at small and at large J2 separated by 
an intermediate quantum paramagnetic phase without magnetic LRO in the parameter 
region J^' < J2 < J2' where J^' ^ 0.4Ji and J^^ a; O.6J1. The GS at low J2 < J^' 
exhibits Neel-ordered magnetic LRO (with a wave vector Q = (vr, tt)), whereas the GS 
at large J2 > exhibits coUinear stripe-ordered magnetic LRO (with a wave vector 
Q = (7r,0))orQ = (0,7r)). 

Given the key role played by quantum fluctuations in determining the gs structure 
of frustrated magnets, it is clearly of central interest to focus special attention on the 
various means by which we may vary or "tune" them. Clearly, as we have already noted, 
an increase in the spin quantum number s is expected to decrease their strength. Thus, 
for example, for the simple case of the isotropic Heisenberg model on the square lattice 
with NN bonds all of the same strength, whereas the quantum fluctuations reduce 
the perfect Neel ordering in the classical case (i.e., s 00) so that the staggered 
magnetization is only about 61% of its classical value for the s = | case as noted 
above, the corresponding reduction in the s = 1 case is less, namely to about 80% of 
the classical value (and see [47] and references cited therein). One of the goals of the 
present paper is to investigate similarly the effect of increasing s for the archetypal J1-J2 
model on the 2D square lattice. In order to do so it is convenient to consider at the 
same time any other means to "tune" the quantum fluctuations. In particular, we note 
that besides changing s or the dimensionality and lattice type of the system, and apart 
from varying the relative strengths of the competing exchange interactions, another key 
mechanism to tunc the quantum fluctuations is the introduction of anisotropy, either in 
real space [48-53] or in spin space [54-57], into the existing exchange bonds. 

Turning flrst to the case of anisotropy in real (crystal lattice) space, we note 
that Nersesyan and Tsvehk [48] have recently introduced and studied an interesting 
generalization of the pure J1-J2 model for the s = | case in order to investigate the 
effects of spatial anisotropy on the quantum fluctuations in the model. This extended 
model, the so-called Ji-J[-J2 model, has been further studied by other groups for both 
the s = I [49-52] and the s = 1 [53] cases. This generalization of the 2D J1-J2 rnodel 
introduces a spatial anisotropy on the square lattice by allowing the NN bonds to have 
different strengths Ji and J[ in the two orthogonal spatial lattice dimensions, while 
keeping all of the NNN bonds across the diagonals to have the same strength J2. In 
previous work of our own [52, 53] on this Ji-J(-J2 model we studied the effect of the 
couphng J[ on the semiclassical Neel-ordered and stripe-ordered phases. For the s — ^ 
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case, we found that the quantum critical points for both of these phases with LRO 
increase as the couphng ratio J[/Ji is increased, and an intermediate phase with no 
magnetic LRO only emerges when J[/Ji ^ 0.6, with strong indications of a quantum 
triple point at J[/ Ji ~ 0.60, J2/J1 ~ 0.33. For J'l/Ji = 1, the results agree with the 
previously known results of the J1-J2 model described above. 

By contrast, for the s — 1 case, we found no evidence for an intermediate phase 
between the Neel and stripe states, as compared with all previous results for the 
corresponding s — ^ case. However, for the s = 1 case we found instead strong evidence 
for a quantum tricritical point at J(/Ji ~ 0.66, J2/J1 ~ 0.35, where a line of second- 
order phase transitions between the quasiclassical Neel-ordered and stripe-ordered phase 
(for J'l/ Ji < 0.66) meets a line of first-order phase transitions between the same two 
states (for J[/Ji > 0.66). For J[/ Ji = 1 the results obviously reproduce those of the 
usual spin-1 J1-J2 model, for which J|/Ji ~ 0.55 ± 0.01. 

Finally, we turn to the main subject of interest in this paper, namely to further 
the study of the 2D spin-1 Ji J2 model on the square lattice by introducing anisotropy 
in spin space. While the influence of the spin anisotropy on the s — j J1-J2 model on 
the square lattice has been studied by various groups [54-57], including ourselves [1], 
relatively little is known for the s — 1 case. 

Our aim here is to further the study of the J^^^-J^^^ model for the s = 1 case, by 
making use of the coupled cluster method (CCM) carried out to high orders by making 
use of supercomputing resources. The CCM (see [2-4] and references cited therein) is one 
of the most powerful and most universally applicable of all known ab initio techniques 
of modern microscopic quantum many-body theory. It is also one of the most accurate 
methods available at attainable levels of computational implementation. We note, in 
the present context, that the CCM is a particularly effective tool for studying highly 
frustrated quantum magnets, where such other numerical methods as the quantum 
Monte Carlo method and the exact diagonalization method are often severely limited in 
practice, e.g., by the "minus-sign problem" for the former case, and the very small sizes 
of the spin systems that can be handled in practice with available computing resources 
for the latter. This is especially true for spin systems with spin quantum number s > |, 
as are of interest here. The CCM has been applied successfully on many previous 
occasions to calculate the ground-state and excited-state properties of a diverse array 
of quantum spin systems [1,4,12,47,52,53,57-71]. 

2. The model 

Exactly as for the -s = | case that we studied earlier [1], the s = 1 J1-J2 Heisenberg 
model employed here has two kinds of exchange bonds, namely the NN Ji bonds along 
both the row and the column directions of the square lattice, and the NNN J2 bonds 
along the diagonals of the squares. The model is then generalized by including an 
anisotropy in spin space in both types of bonds. The anisotropy parameter A is assumed 
to be the same in both exchange terms, thus producing the so-called J^-^^-J^^^ model. 
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with a Hamiltonian given by 

J, + 4 + A<4) 

where the sums over and {{i,k)) run over all NN and NNN pairs respectively, 

counting each bond once and once only. Both exchange couplings arc assumed to be 
antiferromagnetic here (i.e., Ji > and J2 > 0), and henceforth the energy scale is set 
by putting Ji = 1. We shall also only be concerned here with the case A > 0. 

The model has two types of classical ground state (GS), namely a 2;-aligned state 
for A > 1 and an xy-planar- aligned state for < A < 1. Since all directions in the 
xy-plane in spin space are equivalent, we may choose the direction arbitrarily for the 
xy-planar- aligned state to be the x-direction, say. Both of these z-aligned and x-aligned 
ground states further divide into a Neel (tt, tt) state and coUinear stripe states (columnar 
stripe (tt, 0) and row stripe (0,7r)). There is clearly a symmetry under the interchange 
of rows and columns, and hence we only consider the columnar stripe state. The Neel 
states are the classical GS for J2 < | Ji, and the collinear stripe states arc the classical 
GS for J2 > |Ji. The (first-order) classical phase transition between these states of 
perfect classical LRO occurs precisely at J| = |Ji, . 



3. The coupled cluster method 

We now briefly describe the CCM formalism. For further details interested readers are 
referred, for example, to [2-4] and references cited therein. 

In order to use the CCM the first step is always the choice of a normalized model (or 
reference) state |$) which is required to act as a cyclic vector (or, more physically, as a 
generalized vacuum state) with respect to a complete set of mutually commuting multi- 
configurational creation operators, = (Cf)^ that need to be chosen simultaneously. 
The index I here is a set-index that gives a complete labelling of the many-particle 
configuration created in the state C/|$). The requirements on {|$);C/} are that any 
many-particle state can be exactly decomposed as a unique linear combination of the 
states {C/|$)}, together with the conditions, 

{^\C+^O^CY\^) C+ = l, (2) 

[C/,Cj-] = 0=[Cf,C7]. (3) 

The exact many-body gs ket and bra states, whose solutions we seek via the CCM 
calculation at hand, satisfy the respective Schrodinger equations, 

H\^) = E\^), (4a) 
{^\H = E{^\ , (46) 

respectively, with the normalization defined by (^I'l^I') = 1 [i.e., with (^| = 
((^|^'))^^(\1'|], and with |^') itself satisfying the intermediate normalization condition 
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($|\E') = 1 = ($1$). In terms of the set {|$); C*/}, the CCM now employs an exponential 
parametrization for the exact gs ket energy eigenstate, 

that lies at the heart of the method. Its counterpart for the exact gs bra energy eigenstate 
is chosen as 

(^| = ($|^e^, S=l + Y,^iC7. (56) 

The gs CCM correlation operators, S and 5*, contain the real c-number correlation 
coefficients, Si and Si, that need to be calculated. Clearly, once they are known, all 
other gs properties of the many-body system can be derived from them. In order to 
find them we simply insert the parametrizations (l5aj) and (l56j) into the Schrodinger 
equations (Hoj) and (l4B . and then project onto the complete sets of states ($|C7 and 
C/|$), respectively. Completely equivalently, we may simply demand that the gs energy 
expectation value, H = (\E'|if |\E'), is minimized with respect to the entire set {Si, Si}. 
In either case we are easily led to the equations 

($|C7e-^i/e^|$) =0; V/ 7^ , (6a) 
{^\Se-^[H, C+]e^|$) = ; V/ ^ , (66) 

which are first derived using computer algebra and then solved for the set {Si, Si} 
within specific truncation schemes described below, by making use of parallel computing 
routines [72]. Equation (j6aj) also shows that the gs energy at the stationary point has 
the simple form 

E = E{{Si}) = (<l>|e-^i/e^|$) . (7) 

It is important to realize that this bi-variational formulation does not necessarily lead 
to an upper bound for E when the summations for S and S in (ISolb) are truncated, due 
to the lack of manifest Hermiticity when such approximations are made. Nonetheless, 
one can prove [3] that the important Hellmann-Feynman theorem is preserved in all 
such approximations. 

Equation ( l6al) represents a coupled set of nonlinear multinomial equations for 
the c-number correlation coefficients {Si}. The nested commutator expansion of the 
similarity-transformed Hamiltonian, 

e-'He' = H+[H,S] + ^[[H, S],S] + ---, (8) 

and the fact that all of the individual components of S in the expansion of flFoj) commute 
with one another by construction, as in ([3]), together imply that each element of 5* in 
(!5aj) is linked directly to the Hamiltonian in each of the terms in ([8]) . Each of the coupled 
equations (1^ is hence of Goldstone linked- cluster type, which thereby guarantees that 
all extensive variables, such as the energy, scale linearly with particle number, A^. Thus, 
at any level of approximation obtained by truncation in the summations on the index 
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I in fl5 gj) and (l56l) . we may always work safely from the outset in the limit ^ 00 of 
an infinite system, as we do in all our calculations below. It is also important to note 
that each of the linked-cluster equations ( l6a|) is actually of finite length when expanded, 
since the otherwise infinite series of (|8]) will always terminate at a finite order, provided 
only that each term in the Hamiltonian, H, contains a finite number of single-particle 
destruction operators defined with respect to the reference (vacuum) state |$), as in the 
case of our Hamiltonian ([1]). 

We turn now to the implementation of the CCM for quantum spin systems, for 
which it is usually convenient to take the classical ground states as our (initial) choices 
for the model state |$). Hence, we may choose here either a Neel state or a coUinear 
(columnar) stripe state for |$). Each of these can be further sub-divided into a z-aligned 
choice or an xy-plemar (say, x-aligned) choice, which we expect to be appropriate for 
the regions A > 1 and < A < 1 respectively on purely classical grounds. We present 
results in section H] based on all four of these classical ground states as choices for |$). 
In order to implement the CCM computationally it is very convenient to treat the spins 
on every lattice site in any chosen model state 1$) as equivalent. In order to do so 
we introduce a different local quantization axis and a correspondingly different set of 
spin coordinates on each site, so that all spins, whatever their original orientations in 
1$) in the global spin-coordinate system, align along the negative z-direction, say, in 
these local spin coordinates. This can always be done by defining a suitable rotation 
in spin space of the global spin coordinates at each lattice site. Such rotations are 
canonical transformations that leave the spin commutation relations unchanged. In 
these local spin axes where the configuration indices I simply become a set of lattice site 
indices, / {ki, k2, ■ ■ ■ k^}, the generalized multi-configurational creation operators 
Cf are simple products of single spin-raising operators, -^fci -^fe, ' ' ' '^fcm ' "^^ere 

= sl±isl , and (s^, s^, si) are the usual SU(2) spin operators on lattice site k. For the 
quasiclassical magnetically-ordered states that we calculate here, the order parameter 
is the sublattice magnetization, M, which is given within our local spin coordinates 
defined above as 



The CCM formalism is clearly exact if one includes all spin configurations / in the 
expansions (15 a\) and fISfej) of the 5* and 5* operators respectively. However, truncations 
are necessary in practice. Based on a great deal of previous experience, we usually 
employ the so-called LSUBn approximation scheme for s = 1/2 quantum spin systems 
(see [52] and references cited therein), and its so-called SUBn-m counterpart for s = 1 
systems (see [53] and references cited therein). The LSUBn scheme is defined such that 
all possible multi-spin-flip correlations over different locales on the lattice defined by n 
or fewer contiguous lattice sites are retained at the nth level of approximation. For the 
case of spins with s = |, the multi-configurational creation operators, Cf can contain no 
more than one spin-raising operator sj for each lattice site j. However, the number of 




(9) 



k=l 
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fundamental LSUBn configurations for s = 1 becomes appreciably higher than for s = i, 
since each spin on each site j can now be flipped twice by the spin-raising operators, 
so that in this case the multi-configurational creation operators, C/ can contain up to 
two spin-raising operator for each lattice site j. Thus, for systems with s > | it is 
more practical to use the SUBn-m scheme, in which all correlations involving no more 
than n spin flips spanning a range of no more than m adjacent lattice sites are retained. 
Clearly, for spins with s = 1, the SUB2n-n scheme is fully equivalent to the LSUBn 
scheme. More generally for spins with arbitrary spin quantum number s, SUB2s?T,-n 
= LSUBn. In order to keep the number of fundamental configurations from growing 
too quickly with increasing level of approximation we set m = n, and thus we have the 
SUBn-n scheme. The approximation clearly becomes exact as n ^ 00. 

We note that, in general terms, both the LSUBn and SUBn-m truncation schemes 
are systematic localized approximation hierarchies in which the truncation indices 
are physically related to the size of the clusters of spins on the lattice for which 
the multi-spin correlations are explicitly included. Their physical motivation (and 
eventual justification) thus stems ultimately from the localized short-range nature of the 
underlying Hamiltonian (which, in the present case, involves just two-spin interactions 
at NN and NNN distances apart only). The maximum number of spins correlated in 
such clusters is n in both cases. By contrast, the SUBn scheme (which is formally 
equivalent to the SUBn-m scheme in the limit m —>■ 00) explicitly correlates all clusters 
of spins involving no more than n spin-flips, regardless of the spatial separations of the 
spins within the correlated clusters. It is important to note however that in all CCM 
approximations (including the LSUBn and SUBn-m schemes) each correlated cluster 
configuration retained within the correlation operator S of ( !5aj) is actually counted an 
arbitrarily large number of times due to the exponentiated form in which the operator S 
appears in the parametrization fl5aj) . It is precisely the exponential form that guarantees 
the proper counting of arbitrary multiples, at different positions on the lattice, of each 
configuration (and all products of such multiples for different configurations) retained 
in S, considered as independent excitations. Thus, even though, for example, the 
LSUBn and SUBn-m truncation schemes are motivated by the inclusion of the explicit 
correlations within localized clusters of spins only up to a given size, every approximation 
includes configurations in which an arbitrary number of spins (up to all iV — >■ 00 spins) 
are correlated, albeit as (properly counted) products of independent sub-clusters up to 
a given finite size. 

Table [1] shows the number of fundamental SUBn-n configurations for the z- 
aligned and planar x-aligned states in the Neel and striped phases. We see that the 
number of fundamental configurations for the planar model state at the SUB8-8 level of 
approximation is 61904 for the Neel phase and 123471 for the stripe phase. The intensive 
calculations required at even this very high order of approximation are easily practicable 
with relatively modest supercomputing resources. Thus, for example, we employed 200 
processors simultaneously to execute the SUB8-8 calculations using the planar x-aligned 
coUinear stripe state as model state, and with this number of processors it took about 
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Table 1. Numbers of fundamental configurations (U f.c.) retained in the CCM SUBn- 
n approximation for the 2-aligned states and the planar x-aligned states of the s — 1 
jxxz_jxxz ]-QQ(jgi Qj]^ square lattice. 





2;-aligned states 


planar 


x-aligned states 


Scheme 


tl 


f.c. 




tlf.c. 


Neel 


stripe 


Neel 


stripe 


SUB2-2 


1 


1 


2 


3 


SUB4-4 


15 


21 


31 


57 


SUB6-6 


375 


585 


1085 


2131 


SUB8-8 


17864 


29411 


61904 


123471 



six hours to solve the CCM equations fl6aj) and fl6fej) at this level of approximation for 
each value of the anisotropy parameter A in the Hamiltonian ([1]) . 

Clearly, the last step in our calculations is to extrapolate the approximate SUBn-n 
results to the exact, n — >■ 00, limit. We use here for the extrapolations of the raw 
SUBn-n data the same well-tested scaling laws as we used previously in our studies of 
the Ji-J[-J2 model for both the s = | case [52] and the s = 1 case [53], as well as for 
the s = ^ version of the present model [1]. Thus, the scaling law used for the gs energy 
per spin, E/N, is 

E/N = ao + ain'"^ + a2n~'^ , (10) 
and that for the staggered magnetization, M, is 

M = 60 + (&i + M"') . (11) 

In order to have a robust and stable fit to any fitting formula that contains m 
unknown parameters, it is well known that it is desirable to have at least (m + 1) 
data points (the so-called m + 1 rule). Both of our scaling laws (fTOj) and (ITT]) contain 
m = 3 unknown parameters to be determined, and in all cases we have SUBri-n data 
sets with n = {2,4,6,8}. In all our results presented below the SUBn-n results 
are extrapolated to the limit n —>■ 00 using the sets with n = {2,4,6,8} for both 
the 2;- aligned and planar x-aligned states. However, we have also extrapolated E/N 
and M using the sets n = {4,6,8} and n = {2,4,6}. In all cases they lead to 
very similar results, thereby adding credence to their validity and stability. We also 
note that for the corresponding s = 1/2 model we could perform LSUBn = SUBn- 
n approximation calculations for n = {2,4,6,8,10}. This enabled us to perform 
extrapolations using the sets n = {2,4,6,8} and n = {2,4,6,8,10} as well as the 
preferred set n = {4,6,8,10}. Gratifyingly, all sets yielded very similar extrapolated 
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Figure 1. (Colour online) Extrapolated CCM SUBn-n results using the z-aligned and 
planar x-aligncd states for the gs energy, E/N, for the Neel and stripe phases of the 
s — 1 J^^^ ~J^^^ model. The SUBn-n results are extrapolated to the limit 71 cx) 
using the sets n = {2,4,6,8} for both the z-aligned and planar x-aligned states. The 
NN exchange coupling Ji — I. The meaning of the i?max points shown is described in 
the text. 



results, even near phase boundaries and the quantum triple point, which gives us great 
confidence in the accuracy and robustness of our extrapolation scheme. 

4. Results 

Figure [1] shows the extrapolated CCM results for the gs energy per spin, E/N, as a 
function of J2 for various values of A, using both the 2;-aligned and planar x-aligned 
model states. For each value of A two curves are shown, one (for smaller values of J2) 
using the Neel state, and the other (for larger values of J2) using the stripe state as 
CCM model state. As has been discussed in detail elsewhere [3,4,63], the coupled sets 
of LSUBn equations (16 a\) have natural termination points (at least for values n > 2) for 
some critical value of a control parameter (here the anisotropy. A), beyond which no real 
solutions to the equations exist. Thus, for each set of calculations based on one of the 
four CCM model states used, the i^max points shown in figure [U are either those natural 
termination points described above for the highest (SUB8-8) level of approximation we 
have implemented, or the points where the gs energy becomes a maximum should the 
latter occur first (i.e., as one approaches the termination point). The advantage of this 
usage of the -Emax points is that we do not then display gs energy data in any appreciable 
regimes where SUBn-ra calculations with very large values of n (higher than can feasibly 
be implemented) would not have solutions, because of having terminated already. 

All of the curves such as those shown in figure [T] illustrate very clearly that the 
corresponding pairs of gs energy curves (for the same values of A) for the Neel and 



stripe phases cross one another, for both the 2;-aligned (figure 1(a) for all values A > 1) 



and the x-aligned (figure 1(b) for all values < A < 1) cases. The crossings occur 



with a clear discontinuity in slope, which is completely characteristic of a first-order 
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Figure 2. (Colour online) Extrapolated CCM SUBn~n results using the ^-aligned 
and planar x-aligned states for the gs staggered magnetization, M, for the Neel and 
stripe phases of the s = 1 J^^^~J^^^ model. The SUBn-n results are extrapolated 
to the limit n — > cx) using the sets n = {2, 4, 6, 8} for both the z-aligned state and the 
planar x-aligned states. The NN exchange coupling Ji = 1. 



phase transition, exactly as observed in the classical (i.e., s oo) case. Unlike in the 
s = ^ version of this model that we studied earlier [1], there is no indication at all in 
the present s = 1 case of any intermediate paramagnetic phase emerging for any values 
of the parameters J2 and A. Furthermore, the direct first-order phase transition, so 
indicated by our results for the gs energy, between the quasiclassical Neel-ordered and 
coUinear stripe-ordered phases, in both the 2;-aligned and planar x-aligned cases, occurs 
for all values of A > very close to the classical phase boundary J2 = \, the point of 
maximum (classical) frustration. 

We show in figure [2] corresponding indicative sets of CCM results, based on the same 
four model states, for the gs order parameter (viz., the staggered magnetization), to 
those shown in figure [1] for the gs energy. The staggered magnetization data completely 
reinforce the phase structure of the model as deduced above from the gs energy data. 

Thus, let us now denote by Mc the quantum phase transition point deduced from 
curves such as those shown in figure O where Mc is generically defined to be either (a) 
the point where corresponding pairs of CCM staggered magnetization curves (for the 
same value of A), based on the Neel and stripe model states, intersect one another if 
they do so at a physical value M > 0; or (b) if they do not so intersect at a value 
M > 0, the two points where the corresponding values of the staggered magnetization 
go to zero. Clearly, in this generic scenario, case (a) corresponds to a direct phase 
transition between the Neel and stripe phases, which will generally be first-order if the 
intersection point has a value M 7^ (and, only exceptionally, second-order, if the 
crossing occurs exactly at M = 0). On the other hand, case (b) corresponds to the 
situation where the points where the LRO vanishes for both quasiclassical (i.e., Neel- 
ordered and stripe-ordered) phases are indicative of a phase transition from each of 
these phases to some intermediate magnetically-disordered phase. A detailed discussion 



The quantum Jf -J2 spin-1 Heisenberg model on the square lattice 



13 



0.56 
0.55 
0.54 
0.53 
0.52 
0.51 
0.5 
0.49 



: z-aligned □ 

: planar g 

Emeet : z-aligned a 

Emeet ■ P^^nar A 



6 



A- 
/O 
A/ 



'A 



(3 



stripe 



-o 



Neel 



0.5 



1.5 



2 
A 



2.5 



3.5 



Figure 3. (Colour online) Extrapolated CCM SUBn-n results using the z-aligned and 
planar x-aligned states for the ground-state phase diagram of the s = 1 J^-^^-J^^^ 
anisotropic Heisenberg model on the square lattice, for the NN exchange coupling 
Ji — 1. The SUBn-n results for the energy per spin and the staggered magnetization 
are extrapolated to the limit n 00 using the sets n — {2,4,6,8} for both the z- 
aligned and planar x-aligned model states. Ale = magnetization critical point, defined 
in the text. -Emcct denotes the crossing point of the CCM energy curves for the same 
value of A based on the Neel-ordered and coUinear stripe-ordered model states. 



of this order parameter criterion for a phase transition and its relation to the stricter 
energy crossing criterion has been given elsewhere [69]. 



It is clear from figures 2(a) and 2(b) that case (b) above never occurs for the present 
spin-1 model for any values of the anisotropy parameter A or for any values of the NNN 
exchange coupling J2, unlike in the s = ^ version of this model that we studied earlier [1]. 

By putting together data of the sort shown in figures [T] and [2] we can now deduce 
the gs phase diagram of our system from our CCM calculations based on the four model 
states with quasiclassical antiferromagnetic LRO that we have employed. Figure [3] 
shows the zero-temperature gs phase diagram of the 2D s = 1 J^^^-J^^^ model on 
the square lattice for the 2;-aligned and planar x-aligned states, as obtained from our 
extrapolated results for both the gs energy and the gs order parameter. The completely 
independent results from both the energy criterion and the order parameter criterion for 
the phase transition give extremely similar positions for the phase boundary, as one can 
observe from figure [31 Note that the results from using the order parameter criterion 
become increasingly inaccurate for large values of A, and this is why we show them in 
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figure El only out to A < 2. The reason for tliis is simple. Thus, as A — oo, the order 
parameters M —>■ 1 for both the Neel-ordered and coUinear stripe-ordered phases, and it 
becomes increasingly difficult to determine the point where they cross, since the angle of 



their crossing becomes vanishingly small. This effect can clearly be seen in figure 2(a 
where it has clearly become acute even for values of A as small as about 2. On the 
other hand, the energy criterion correspondingly becomes more accurate as A ^ oo, as 



one may observe from figure 1(a) Thus, figure [3] clearly shows that the phase boundary 
approaches the classical line J| = 0.5 as A — * oo, as expected in this Ising-like limit. 

Our results certainly provide very clear and consistent evidence that there exists 
no intermediate phase. Thus, the curves for the order parameters of the Neel and stripe 
phases always meet at a finite value and the corresponding curves for the gs energies of 
the two phases intersect with a discontinuity in slope, for both the 2;-aligned and planar 
x-aligned states, for all values of the anisotropy parameter A. All of the evidence clearly 
points towards a first-order phase transition between the two phases. 

We note also that the z-aligned and a;|/-planar-ahgned phases meet precisely at the 
isotropic point A = 1, just as in the classical case, and exactly as expected. However, 
this does provide a consistency check on our independent numerical calculations for the 
two phases. The case A = 1 obviously reproduces the usual (isotropic) J1-J2 model. 
Thus, at A = 1, we find J2 = 0.55 ± 0.01 which, very encouragingly, is the same value 
we found [53] for the s = 1 Ji-J[-J2 model in the spatially isotropic limiting case when 
J(/Ji = 1. We also note that in the present spin-1 quantum model, the isotropic point 
A = 1 is precisely the point at which the boundary between the two quasiclassical 
phases deviates most from its classical position at J| = | for all values of A > 0. Our 
calculations also indicate that at the isotropic XY point of the model (i.e., where A = 0) 
the phase boundary is at J| = 0.50 ± 0.01. 



5. Discussion 



Our results have clearly shown in detail how the quantum fluctuations present in the 
spin-1 J1-J2 model on the infinite square lattice are diminished by varying the spin 
anisotropy parameter A away from the Heisenberg isotropic point A = 1 in either 
direction. This is precisely as was observed previously [1] for the spin-i version of the 
same model, and as was to be expected. However, unlike what would be predicted by 
lowest-order (or linear) spin- wave theory (LSWT) [6], for example, we can now conclude 
with confidence from our results that no such intermediate disordered phase as the one 
that we observed in the spin-^ version of this model between the two quantum triple 
points at (A'^ = -0.10 ± 0.15, J^/Ji = 0.505 ± 0.015) and (A'^ = 2.05 ± 0.15, J^/Ji 
= 0.530 ± 0.015), exists for the spin-1 version, for any values of the parameters J2/J1 
and A. 

In the context of a spin-wave theory (SWT) treatment of the isotropic J1-J2 model 
on the square lattice, LSWT predicts that quantum fluctuations can destabilize the 
classical GS with LRO, even at large values of the spin quantum number s, for values 
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of the frustration parameter J2/J1 around 0.5. For the spin-| case the range of values, 
a'^^ < J2/J1 < o!^^ 1 for which a magnetically-disordered phase thereby occurs is predicted 
by LSWT to be given by ~ 0.38 and a^'^ 0.52. These values may be compared 
to our own predictions [1] of a^^ = 0.44 ± 0.01 and a'^^ = 0.59 ± 0.01. For the spin-1 
case LSWT predicts a narrower, but still non- vanishing, strip of disordered intermediate 
phase in a range with a*^^ 0.47 and a'^'^ 0.501, whereas we predict with confidence 
that the disordered phase simply does not exist as a GS in this case. 

The discrepancy between our results and those of LSWT for the spin-1 case are 
undoubtedly due to the shortcomings of LSWT. Thus, while LSWT can work reasonably 
well in the absence of frustration (e.g., for the isotropic J1-J2 model here when J2 = 0, 
that represents the Heisenberg model with only NN interactions), in the presence of 
frustration it consistently overestimates the effects of quantum fluctuations. This effect 
worsens as the frustration (here measured by the ratio J2/J1) increases. 

Thus, Igarashi [73] has shown explicitly for the J1-J2 model by going to higher 
orders in SWT (i.e., by calculating higher-order terms in the 1/s power expansion), 
that while the series seems to converge for values J2/J1 < 0.35, the second-order 
corrections grow so large for values J2/J1 > 0.4 that no prediction based on LSWT, or 
even on higher-order SWT, in this region (e.g., about the appearance of an intermediate 
magnetically-disordered phase near J2/J1 ~ 0.5) should be relied upon. Furthermore, 
he showed that the effects of the higher-order correction terms to LSWT make the 
Neel-ordcrcd state more stable than predicted by LSWT. 

Relatively little attention has been paid by other authors to the (pure, isotropic) Ji- 
J2 model at higher values of the spin quantum number, s > |. We note, however, that 
Cai et al. [74] have also recently postulated the possible existence of an intermediate 
phase between the quasiclassical Neel-ordered and coUinear stripe-ordered phases for 
the spin-1 model. More specifically, they hypothesize an intermediate valence-bond 
solid (VBS) ground state (GS) for the spin-1 isotropic J1-J2 model at or near the point 
of maximal classical frustration where J2/J1 = 0.5. Their evidence is indirect and is 
based on a trial variational state of VBS type, which is an exact GS of a related spin-1 
model Hamiltonian, and on a pscudopotential approach to extend it to the actual spin- 
1 J1-J2 model. They express the dual hopes that this trial state might capture the 
main character of the disordered phase that they thereby predict for the fully frustrated 
case, and that accurate numerical methods, such as those considered here, might verify 
the existence of this postulated intermediate phase. Such variational analyses, based 
on physically motivated trial states, are always of interest, but have a very chequered 
history of success in the field of highly correlated spin- and electron-lattice systems. In 
the present case we stress again that our own detailed numerical analysis provides no 
evidence at all for the existence of such an intermediate magnetically-disordered VBS 
phase as postulated by Cai et al. [74]. 

In the same context, we note too that in earlier work Read and Sachdcv [75] have 
applied a large- expansion technique based on symplectic Sp(A^) symmetry to the 
isotropic J1-J2 model. They found that the method, which can itself be regarded as akin 
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to a 1/s expansion, predicts an intermediate phase (with VBS order) for smaller values 
of s, but that this phase disappears for larger values of s where they predict instead a 
first-order transition between the Neel and stripe phases. All of these qualitative results 
for the pure J1-J2 niodel are in accord with our quantitative predictions. 

We note that the results presented here for the spin-anisotropic spin-1 Jf^^-J^^^ 
model are also fully consistent with our own previous results [53] for the spatially- 
anisotropic spin-1 Ji-J[-J2 model discussed in section [1] above, for which we also found 
no evidence for an intermediate disordered phase between the quasiclassical Neel and 
collinear stripe phases with LRO. However, whereas for the spin-1 Ji-J[-J2 model we 
found strong evidence for a quantum tricritical point at {J[/ Ji ~ 0.66, J2/J1 ~ 0.35) 
where a line of second-order phase transitions between the Neel-ordered and the collinear 
stripe-ordered states (for J[/Ji < 0.66) meets a line of first-ordered phase transitions 
between the same two states (for J[/Ji > 0.66), we find for the present spin-1 J^^^- 
jxxz ]2iodel that the phase transition between these two states is first-order for all 
values A > 0. Clearly, these two sets of results are in complete agreement with one 
another at their common point of overlap, when J[ = Ji and A = 1. 

At the XY isotropic point (A = 0) of the present spin-1 J^^^-J^^^ model we 
predict that the phase boundary occurs at a value J|(0) = 0.50±0.01. It is interesting to 
note that our previous results for the spin-1 version of the model [1] showed a quantum 
triple point (QTP) at (A'^ = -0.10 ±0.015, = 0.505 ± 0.015). Clearly our results for 
this spin- 1 case are consistent with this lower QTP occurring exactly at the XY isotropic 
point (A = 0) and also at the point of maximum classical frustration, J2 = |- Similarly, 
in the present spin-1 case our results are consistent with the phase boundary at the XY 
isotropic point also occurring at the point J2 = |- It would seem likely, therefore, that 
for both the cases of spin-i and spin-1 particles the corresponding quantum J{^^-J^^ 
model has a special behaviour at the point J2/J1 = \ where the classical frustration 
is greatest. Our results indicate that a more detailed investigation of this case might, 
therefore, be worth undertaking for general values of the spin quantum number s. 

Although there is very little other accurate numerical work for the present 
model against which to make comparisons, there have been several previous detailed 
comparisons, for example, of CCM results with those from the exact diagonalization 
(ED) of finite spin-lattices for some particular models. One such example [65] is the 
spin-1 J- J' (or zigzag) model on the square lattice which contains two kinds of NN 
isotropic Heisenberg interactions, of strength J and J' respectively, such that each 
square plaquette contains three J-bonds and one J'-bond, with the J'-bonds arranged 
in a regular zigzag fashion such that every lattice site on the square lattice is joined 
to only one J'-bond. An alternative but equivalent description of the model is that it 
interpolates between a honeycomb and a square lattice, such that the J-bonds join NN 
lattice sites on the honeycomb lattice, and the J'-bonds join sites across only one of the 
main diagonals of each hexagon, such that when J = J' the model is equivalent to the 
NN isotropic Heisenberg model on the square lattice. 

ED calculations were performed for the above model [65] for lattices with up to 
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N = 32 sites. In general terms it was found that the CCM results for the model at 
attainable levels of implementation (viz., using the LSUBn approximation with n < 8 
agree well with the extrapolated (A^ 00) ED data. The CCM is particularly good, 
however, at describing both the dimerized and the helical gs phases that this system 
can support. For the latter phase the ED results lie appreciably above those from the 
CCM. This is because the energies for the small lattices able to be considered do not fit 
well to the known theoretical finite-size scaling law in this regime. It is no surprise that 
finite-size effects for systems with an incommensurate helical spin structure are larger 
than for systems with Neel order or that are ordered with dimerized spin pairs. 

Similar conclusions were also drawn for comparions of CCM and ED results for 
extensions of the above spin-i J- J' model to both (a) the anisotropic Jxxz^J'xxz 
model [57] where both bonds contain an Ising anisotropy of precisely the sort considered 
in the present paper; and (b) the case where the spin quantum number s > | [68]. For 
the latter case of the spin-1 J- J' model, calculations were performed using both the 
CCM in the SUBn-n scheme with n < 6 and the ED technique on lattices of sizes 
N < 20. Again, the resulting finite-size ED extrapolations remained quite poor, and 
only allowed some qualitative conclusions to be drawn, whereas results from the CCM 
were seen to be much more robust and more rehable. In no case, however, did the CCM 
and ED results conflict with each other. 

Another model where ED and CCM results have been compared is the pure 
(isotropic) spin-| J1-J2 model on the square lattice [76]. Again, the ED results (with 
< 32) were found to provide a good qualitative check of the CCM data for LSUBn 
calculations performed with n < 8. Finally, for the spin-i version of the present 
anisotropic Jf^^'^- J^^^ model, we [77] have also compared the CCM results with those 
from ED calculations on finite-sized lattices of size A?" = 36 = 6 x 6 sites (with periodic 
boundary conditions imposed) . In this case too the ED data are best used to complement 
the CCM results. On the basis of all the above evidence we expect that the same will 
hold true for the spin-1 version of the model studied here. Since the number of basis 
states increases roughly as 3^ for the spin-1 case, by comparison with 2^ for the spin-| 
case, ED calculations for the present model would be limited to lattices of sizes = 16 
and A^ = 20. The next biggest lattice that preserves the full lattice symmetry has 
N = 26 sites in this case, and an ED calculation of this size for the spin-1 model is 
probably beyond the limits of presently available computing power. With only such 
limited data the ED finite-size extraploation would again be bound to remain poor, as 
seen in the previous work cited above, and we fully expect that the CCM results would 
again prevail even if ED results were available for the present model. 

Finally, we note that our analysis and conclusions have relied heavily on two of 
the unique strengths of the CCM, namely its ability to deal with highly frustrated 
systems as easily as unfrustrated ones, and its use from the outset of infinite lattices. 
These, in turn, lead to its ability to yield accurate predictions for the locations of phase 
boundaries. Our own results for the gs energy and staggered magnetization provide a 
set of independent checks that lead us to believe that we now have a self-consistent 
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and coherent description of these challenging anisotropic and frustrated J^-^^-J^^^ 
systems for both the spin-i and spin-1 cases. 
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